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Abstract. We study a planar two-temperature diffusion of a Brownian particle in 
a parabolic potential. The diffusion process is defined in terms of two Langevin 
equations with two different effective temperatures in the X and the Y directions. 
In the stationary regime the system is described by a non-trivial particle position 
distribution P(x,y), which we determine explicitly. We show that this distribution 
corresponds to a non-equilibrium stationary state, characterised by the presence of 
space-dependent particle currents which exhibit a non-zero rotor. Theoretical results 
are confirmed by the numerical simulations. 
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1. Introduction 

The idea of physical systems characterised by two different temperatures has been 
proposed a long time ago for the models of spin-glasses and neural networks with 
partially annealed disorder [H |2l |3l H]. In these models, the two temperatures T\ 
and T2 are related to two different degrees of freedom, which are evolving at two 
essentially different time scales. As an example, one may consider a system in which 
the fast spin variables are connected with the thermal bath kept at the temperature 
Ti, while the slow spin-spin coupling parameters are connected with another thermal 
bath maintained at the temperature T 2 . It can be easily shown that in the stationary 
(non-equilibrium) state the statistical properties of such systems are described by the 
usual replica theory of disordered systems with a finite value of the replica parameter 
n = T1/T2 (see also [5]). Unfortunately, generalisation of this idea to the case when 
dynamics of two types of degrees of freedom is characterised by two comparable (or 
equal) time scales turned out to be rather problematic: it seems that there is no 
generic explicit expression for the stationary probability distribution function which 
would generalise the Gibbs distribution of the equilibrium case T\ = T 2 [6]. However, 
there is a particular case for which one can find an explicit and a rather non-trivial 
expression for the stationary distribution function. Namely, this is the case when the two 
dc grees of freedom x and y related to the thermal baths with the temperatures T x ^ T y 
respectively, experience a potential which is a quadratic function of x and 2/ [HJ [T] . During 
last decade theoretical investigations of such type of systems were mostly concentrated 
on the studies of nonequilibrium fluctuations and energy transfer [8]. Recently this type 
of model was studied both theoretically j9] and experimentally [10] from the point of 
view of the entropy production and memory effects. In this paper, keeping in mind 
putative experimental realisation of such a type of systems, we are going to discuss the 
two-temperatures situation reformulated in terms of the two-dimensional diffusion of a 
Brownian particle in a parabolic potential. The diffusion process is defined in terms 
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of Langevin dynamics with two different effective temperatures in the X and the Y 
directions. In the stationary state this system is described by a non-trivial distribution 
function P(x,y), which can be computed explicitly. Unlike for the equilibrium case 
(T x = T y ), this non-equilibrium stationary state is characterised by the presence of 
nontrivial space dependent particle's flows y). Moreover, these flows exhibit a 
"symmetry breaking" rotor, S(x, y) = V x y) (directed perpendicular to the (X, Y)- 
plain), the sign (or the direction) of which is determined by the temperature difference 
(T x -T y ). 

The paper is organised as follows. In Section |2] we define our model and present the 
explicit solution for the stationary particle's probability distribution function P(x,y). 
In Section [3] we compute putative "observable" quantities of the system, such as the 
variances of the particle displacements in the X and the Y directions, the rotor S(x, y) 
of the particle's flows as well as the average rotation velocity. In Section H] we report the 
results of the numerical simulations and compare them with our analytical predictions. 
Finally, in Section [5] we conclude with a brief recapitulation of our results. 



2. The model 



We consider stochastic, over-damped Langevin dynamics of a particle moving in a 
two-dimensional space in a presence of an external potential U(x,y). The particle 
instantaneous position p{t) is defined by projections on the X and the Y axes, x(t) 
and y(t), respectively. The time evolution of x(t) and y(t) is described by following 
equations: 

± X (t)=-^U(x,y) + Ut) 

(1) 

d , s d 
dt m = ~dy- 

Here £ x ,y(t) is anisotropic stochastic noise, with zero mean and correlation function 

(UW)> = 2T a 5 a> p6(t-t') , (a, (3 = x, y) (2) 

where T x and T y are two different" temperatures" and U(x, y) has the following parabolic 
form : 

U{x,y) = ^x 2 + ^y 2 + uxy, (3) 

The shape of the potential is controlled by the parameter u. To keep the particle 
localised near the origin, we have to impose the constraint |u| < 1. This follows from 
the requirement that both eigenvalues of the potential, Ai^ = I ± n, must be positive; 
in the case \u\ > 1, there is a direction in the plane (x, y) at which the potential U(x, y) 
has a negative curvature which allows the particle to escape to infinity. 
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In the stationary regime, the probability distribution function P(x, y) of the particle 
position obeys the stationary Fokker-Planck equation : 







T x dP ^ y) + P(x,y) dU( ?' y) ] + l\T y dP ^ V) + P(z,y) dU< ?> y) = (4) 



OX L OX OX 



dy 1 y dy ' dy 



In the trivial isotropic case, T x = T y = T, the solution of the above equation is simply 
the equilibrium Gibbs distribution P iso (x,y) oc exp{ — ^ U(x,y)}. 

One can easily show that in the generic anisotropic case with arbitrary T x and T y , 
the solution of the stationary equation fll]) reads: 

P(x,y) = Z _1 exp|--7ix 2 - ^72 - uj 3 xy} (5) 
where the following shortenings have been used 



7l = — 2 v x 2*i (6) 

11 T x T y (l + u 2 A 2 ) ' 1 } 



T x + \u 2 (T x -T y) 

T. + ln 2 ^-^ .) 
T^TJl + u 2 ^ 



79 = y 1 2 v y m 

rp rp I -i i „,2 A 2\ ' 
- a;- 1 1/ 1 



and 



73 ~ 2T a j;(i + w 2 A 2 )' 1 ' 



^y/TyPx 

Further on, Z is the normalisation constant (the "partition function"), defined as 




alx dy exp j - -j t x 2 - \p 2 y 2 - ui3 xy } 



2 T,r„( 1 + „W) 

V 1 — it 

One immediately observes that Z exists, so that the system has the stationary solution, 
only for \u\ < 1. 

3. The observable quantities 

3.1. Variances of particle positions 

Using the above probability distribution function we can straightforwardly calculate the 
variances of the particles position with respect to the X and the Y axes: 
2 = T x + \u\T -T x ) ^ 
1 — u 

(y 2 ) = T y +1 f( T *~ T »\ ( i2) 

1 — u z 
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The characteristic quantity, which can serve as the measure of anisotropy in the system 
under study, is defined as the ratio of these two quantities : 



g{T y /T x ; u) 



(x 2 ) 



2 + u 2 {T y /T x -l) 



(13) 



(y 2 ) 2T y /T x + u 2 (l-T y /T x ) 

In the trivial decoupled case, u = 0, we find g(T y /T x ; 0) = T x /T y , while in the isotropic 
case, T x = T y we have g(l; u) = 1 for all values of the coupling parameter u. Note next 
that in the strongly anisotropic case, e.g., when T y /T x ^> 1, one has 

,2 

-T„ (14) 



(x 2 ) 

(y 2 ) 



u 



2(1 -u 2 ) y 

2-u 2 
2(1 -u 2 ) y 



(15) 



9 



u 



(16) 



In other words, in the strongly anisotropic case the values of both (x 2 ) and (y 2 ) are 
defined by the largest T, while the value of the ratio g = (x 2 ) / (y 2 ) becomes a T- 
independent constant. 



3.2. Mean rotation velocity 

In the stationary case the current j = (j x ,j y ) is defined as follows: 



3x 



^ dP(x,y) .dU(x.y) 
T, ^^ + P(x,y)- 



jy — T t 



X a 

OX 

dP{x, y) 



dy 



P(x,y) 



dx 

dU{x,y) 
dy 



Using eqs.dn])-© we obtain 



1 - T x n/i)x + u(l -T x >y 3 )y P{x,y) 



(17) 
(18) 

(19) 



Jv 



'1 -T y j 2 )y + u(l -T w 7 3 )x P{x,y) 



(20) 



Note that in the isotropic case, T x = T y = T, we have 7i = 72 = 73 = 1/T, so that 
j = 0. In the anisotropic case T x 7^ T y the above non-trivial pattern of currents can be 
characterised in terms of the rotor: 

d d 

S(x,y) = V xj(x,y) = —j y -—j x (21) 

In general, the rotor S(x, y) is a rather complicated function of two variables x and y, 
but it is remarkable that the function S(x, y) has a non-zero (and very simple) value at 
the origin at x = y = 0: 



5(0) 



r-l U T x 



u(T x -T y ) 1?J Z = — 

x y 



y 




(22) 
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Note that this quantity changes sign from minus ("left rotation") at T y > T x , to plus 
("right rotation") at T y < T x . 

Due to the presence of a non-zero particle's current rotor, one finds that the mean 
particle's rotation velocity is also non-zero. Indeed, for a given value of the particle's 
linear velocity v located in the point r on the two-dimensional plane, its angular velocity 
is 

u(t) = 1 (v x r) (23) 

where (v x r) is the vector product directed along the z-axis. Thus, the mean rotation 
velocity (oj) in the limit of an infinite observation time can be defined as follows: 

(to) = lim - / dtu(t) (24) 



Changing averaging over time by averaging over ensemble (which will be justified in 
what follows by numerical simulations) we get: 

//•27T /*00 
d 2 r — (j x r) = J d(j) J dr (j x sin <\> - j y cos 4>) (25) 

Here the average current j is defined in eqs. ( |T9l) - ( 1201) . According to eq. (j5J) , the probability 
distribution function P(r, 0) can be represented as follows 

P(r,0) = Z- 1 exp{-^r 2 M/(0)} (26) 

where 

W{4>) = 7icos 2 (0) + 72sin 2 (0) + m 73 sin(20) (27) 

Substituting the explicit expressions for the components j x and j y of the current, 
eqs. (fT9l) - (l20jl . and using eqs.fjS])-©, we get 

(j x sin0-j,cos0) = M(T ^~ Tx) rW(0) exp{-^ r 2 *(0)} (28) 
Substituting eq.( )28|) into eq.(|25l) and performing simple integrations we obtain 



<"> = " a Vt^ < 29 » 

where the parameter A is defined in eq.Q. 

One can easily prove that the maximal value of the mean angular velocity is 
(u)max = 1) and it is achieved either in the limits A — > — oo (which corresponds to 
T y — > for finite T x ) or in the limit A — > +oo (which corresponds to T x — >■ for a finite 
T y ), and the value of the coupling parameter u = 1/yfK — > 0. 
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4. Numerical simulations: Brownian dynamics 

To verify our analytical predictions and the underlying assumption that the time- 
average can be replaced by the ensemble average, we perform numerical simulations 
of appropriately discretised Langevin equations eqs.flTJ). Substituting the potential 
U(x,y) = \x 2 + \y 2 + uxy into eqs.flTJ) we first write these equations explicitly: 

±(t) = -x-uy + £ x (t) 
y{t) = -y-ux + £ y {t) , 

where the variances of the thermal noise components are defined by = 2T X , 

= 2T y and (Uy) = 0. 

Discretising Eq. (1501 with a time step At, we have : 

x(t + At) = x{t) - At(x + uy) + g x {t)y /TLAt 

y(t + At)=y(t)-At(y + ux) + g v (t)y/2TJM, 1 ' 

where g x (t) and g y (t) are delta-correlated random numbers with Gaussian distribution of 
unit half- width, At <C 1, \V2T x>y At <C 1, which are the conditions of a smooth motion. 
In that case for a free motion of a particle (U(x, y) = 0) which starts at the origin 
((x(0) = 0, y(0) = 0)), the diffusion coefficients are D a = T a , a = x,y and the variances 
of the displacement are given by (x 2 (t)) = 2T x t and (y 2 {t)) = 2T y t. In the case of the 
symmetric potential (u = 0), one has in the stationary regime (x 2 ) = T x and (y 2 ) = T y , 
independently of t. For asymmetric potential m ^ 0, we will compute the mean angular 
velocity (cu) given in eq. (|24j) and the measure of anisotropy g(T y /T x ,u) = (x 2 ) / (y 2 ) 
that is described by eq. (TT3]) . 

The numerical simulation has been done for the time step At = 0.001. The 
averaging has been performed over the total time period r = 10 6 time units, the 
numerical inaccuracy has been evaluated by splitting the whole time interval into 10 
sub- intervals. In Fig. QJa),(b) we plot numerical results for the ratio of variances 
(x 2 ) I (y 2 ) and for the mean angular velocity (a;), calculated as the time-average of 
u(t), as functions of u for T x — 1, T y — 1,2,4,6,8,10. For comparison we also 
show our analytical predictions in eq.( lT3j) and eq.flMJ), respectively, and find a perfect 
agreement. This justifies the replacement of the time-average by the ensemble average 
in our analytical calculations. 

Further on, in Fig. 12(a), (b) we plot the same quantities as functions of T y /T x (with 
T x = 1) for u — 0, 0.2, 0.4, 0.6, 0.8, 0.9. We again observe a very good agreement between 
our numerical and analytical results. 

5. Conclusions 

In the present work we studied a simple stochastic "toy model" with only two 
degrees of freedom which are connected to two thermostats maintained at two different 
temperatures T x and T y , respectively. The model describes diffusion of a particle on a 
two-dimensional plane in a presence of a parabolic potential such that the stochastic 
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U U 



Figure 1. (a) The ratio (x 2 ) / (y 2 ) of variances of particle's displacements along the 
X and the Y-axes vs the parameter u. (Symbols and the color-code is as in (b) ); (b) 
the mean angular velocity (u>) as a function of u for different T y /T x = 1,2,4,6,8, 10. 
Solid lines arc our predictions in eqs.(fT3)) and 




012345 01234 



Ty /T x Ty /T x 

Figure 2. (a) The (x 2 ^ / (y 2 ^ and (b) the mean angular velocity (uS) as functions 
of T v /T x for different values of u = 0, 0.2, 0.4, 0.6, 0.8, 0.9. Solid lines define our 
theoretical predictions in cqs. ([T3| and (|24j) . and the symbols denote the results of 
numerical simulations. 

noises in the X and the Y directions have different strength (T x and T y , respectively). 
We determine the stationary state probability distribution function for the position of 
the particle. Despite its relatively simple structure, it turns out to be rather non-trivial, 
revealing interesting qualitative physical phenomena. In particular, in the stationary 
state one finds a rather sophisticated pattern of particles' density currents (which would 
be identically equal to zero in the equilibrium case) characterised by the non-zero rotor. 
Moreover, due to the presence of this flux rotor one observes the phenomenon which 
could be interpreted as a "spontaneous symmetry breaking", namely one finds non- 
zero value for the average particle's rotation (around the origin) velocity. This value is 
proportional to (T y — T x ), eq.(!29l). being positive (left rotation) for T x < T y and negative 
(right rotation) for T x > T y . 

It should be stressed, however, that except for recently proposed two-temperature 
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electric analog system [10], for the moment the considered model has no experimental 
realization. Thus, the aim of the present work is somewhat provocative: we would like 
argue that the systems of such type are sufficiently interesting to stimulate investigations 
for their "hardware" implementations. We also believe that modifications of our toy 
model towards a system that could be realised in practice and at the same time would 
not loose its interesting behavior (rotation), is possible. 
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